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DATA-REDUCTION TECHNIQUE FOR DETERMINING 
RAPIDLY RISING HEATING RATES TO 
THERMALLY THICK WALLS 

By James L. Dillon, Robert L. Wright, 
and G. Louis Smith 
Langley Research Center 

SUMMARY 

The integral method of determining aerodynamic heat flux from temperatures mea- 
sured at four locations within a thermally thick wall has been analyzed for rapidly rising 
heating rates and large thermal as well as physical thicknesses. A restricted-exponential 
function has been developed for curve fitting the four values of c w (T)-^ through the wall 
(where c w denotes specific heat of wall material; T, temperature; and t, time). The 
restricted -exponential function more accurately describes the variation of c w (T)-^ 
through the wall than the previously used polynomial curve fit and improves the accuracy 
of the integral heating method for determining rapidly rising heating rates to walls of 
large thermal thickness (large temperature gradients through the wall). Results are pre- 
sented for simulated flight measurements for both smooth temperature data and scattered 
temperature data (standard deviation of 15° R, or 8.3° K) for two arrangements of sensor 
locations through the wall. 


INTRODUCTION 

The methods and techniques of analyzing telemetered instrument measurements are 
a vital part of a flight experiment. Free -flight experimental data are, by nature, very 
sensitive to the data-reduction technique utilized. This is due to the test time, sampling 
rates of the measurements, environmental variations throughout the test period, and many 
possible sources of inaccuracy in the telemetered data. 

Numerous methods have been developed for determining experimental heating results 
when test conditions require thick-wall calorimeters. These methods are discussed in 
reference 1. However, only the methods of references 1 to 7 deal with the nonlinear 
heating problem involving temperature -dependent thermal properties of the wall. One of 
these, the integral method as applied in reference 2, has been used successfully in evalu- 
ating the heat flux to a thermally thick (0.508 -cm -thick) heat shield of the Fire spacecraft 


; 
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(refs. 2 and 3). The adaptation in reference 2 of the integral method for integrating 
c w (where c w denotes specific heat of the wall material; T, temperature; and t, 

time) through the wall to determine the heating rate incorporates a third-order polyno- 
mial expression to fit experimental values of c w obtained from temperature mea- 
surements at four locations through the wall. 

In the present paper an assessment of this adaptation for walls of greater physical 
and thermal thickness than the Fire heat shield and a rapidly rising heat flux indicated 
that the accuracy could be improved by incorporating a restricted-exponential function 
to fit the experimental values of c w 4^ through the wall. Data are presented for two 
arrangements of thermocouple locations, with both smooth and randomly scattered tem- 
perature histories. 


SYMBOLS 

A n coefficients in temperature expression (where n = 0, 1, 2, 3, 4, 5, 6) 

a n coefficients in exponential function defining c w (where n = 0, 1, 2, 3) 

c w specific heat of wall material 

k conductivity 

l wall thickness 

q heat flux (heating rate) 

S sum of the squares of the residuals 

T temperature 

t time 

x depth 

p w material density 

Subscript: 


i 


thermocouple location 


ANALYSIS 


Integral Method 

For a nonablating, homogeneous wall with no heat losses from the back surface, the 
one -dimensional (no temperature gradient along the wall) heat flux to the front surface is 
equal to the rate of change of heat stored within the wall. The heat flux can be expressed 
by 


q = p w Iq c w( t )|F ^ (1) 

The integral method for solving equation (1) requires the definition of the gradient 
of c w (T)-^- through the wall and the integration of the area under the curve to produce 
the heat flux. The accuracy in determining heat flux by this method depends on the defi- 
nition of this gradient (analytical curve fit). 

Analytical model. - Evaluation of the integral for a very large temperature gradient 
through the wall was made by analytically simulating a flight heating-rate history and 
comparing the heating rates computed by the integral method with the simulated (known) 
heat input. The heating-rate history of figure 1 (typical of the side -wall heating on a 
slender -cone reentry vehicle at zero angle of attack) was applied in a direct thick-wall 
transient -heating program to a beryllium wall 1.524 cm thick. This program utilizes an 
explicit finite -difference computational technique to determine the temperatures within a 
wall for prescribed heating rates. Thermal measurement locations were assumed to be 
equally spaced within the wall at depths of 0.0254, 0.5080, 1.0160, and 1.5176 cm from 
the heated surface. The thermophysical properties of beryllium shown in figure 2 were 
used in the investigation. A detailed description of the thick-wall heating program may 
be found in the appendix to reference 1. 

The accuracy of the temperatures computed by this program is dependent upon the 
size of the time interval used in the computations and the number of elements into which 
the wall is divided. An accuracy analysis of the thick-wall program for a 1.524 -cm -thick 
wall divided into 11, 30, and 70 elements and for a 0.05-second time interval is presented 
in reference 1. For the heating-rate input of figure 1, a time interval of 0.05 second, and 
a wall division of 30 elements, as used in these computations, the temperatures at each 
location in the wall are estimated to be accurate within ±1 percent. The temperature 
histories for the four assumed measurement locations are shown in figure 1. 

Polynomial curve fit .- The technique (developed in ref. 2) for evaluating the integral 
of equation (1) requires temperature histories at four thermocouple locations in the wall. 
One of these locations is as near to the heated surface as possible, and another at the 
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inside surface; the other two are equally spaced in the wall. The temperature history 
for each location (denoted by subscript i) is fitted with a sixth -order polynomial using 
the Chebyshev technique: 


T i(t) = A 0>i + A 1?i t + . . . + A 6 i t 6 


( 2 ) 


where i = 1, 2, 3, 4. 

For any selected instant of time, the temperature derivatives with respect to time 
■231 are determined from the sixth-order polynomial, and the specific heat of the mate- 
rial c w is obtained as a function of temperature for each thermocouple location. The 
distribution through the wall of the product c w (T)-|~£ (defined by the four locations in the 
wall) is then fitted with a third-order polynomial. The polynomial curve is analytically 
integrated (in closed form) from the heated surface to the inside surface and multiplied 
by the density of the material to yield the rate of change of stored heat in equation (1). 

Temperature histories (fig. 1), simulating measurements from a flight experiment, 
were utilized in the integral method from reference 2 to compute heating rates for com- 
parison with the exact heating-rate input of figure 1. Since no errors of data scatter 
were introduced, the comparison of the simulated and computed heating rate, presented 
in figure 3, provides an assessment of the accuracy of the data-r eduction method. The 
computed heating rates are consistently higher than the simulated heating rates. Although 
the agreement is reasonably good at the lower heating rates, the accuracy of the method 
of reference 2 is noted to deteriorate for very large temperature gradients through the 
wall (higher heating rates). 

The cause of the inaccuracy was investigated and was determined to be the shape 
of the analytical curve used to fit the values of c w through the wall. Figure 4 pre- 
sents the variation through the wall of c w at three intervals of time. The symbols 
indicate the simulated experimental values of c w at the four measurement locations 
which were obtained by differentiating the sixth-order polynomials fitted to the tempera- 
ture histories and multiplying by the appropriate value of c w (T). The symbols are in 
excellent agreement with the solid line which shows the "actual" distribution of c w 
through the skin, as provided by the direct computation of interior wall temperatures for 
the exact heating rate. At the early time (t « 2 seconds), the third-order polynomial 
which is fitted to the four simulated experimental data points provides adequate agree- 
ment. At the later times (t ~ 6 and 7 seconds) when the heating rate is rising rapidly, 
the shape of the third-order -polynomial curve differs significantly from the actual dis- 
tribution, so that the integrated areas under the curves, and hence the heating rates, do 
not agree, as indicated in figure 3. 
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Ex ponential curve fit . - Since the accuracy of heating rates computed by the integral 
method of reference 2 is dependent on the curve fit of the four measured values of c w , 
numerous analytical expressions (polynomials and exponentials of various order) were 
systematically investigated to find a method which would define the distribution of c w 
through the wall more accurately than does a third-order polynomial. Of the expressions 
investigated, an expression of the form 


8T aQ+a-jX+agX^i-agx 3 
Cw ~dt = e 


( 3 ) 


was found to provide the most accurate definition of the gradient of c^ Figure 5 
shows the curve fit of the four simulated points obtained with the exponential function. 
Comparison with figure 4 shows that the exponential -function curve more closely approx- 
imates the shape of the actual curve than does the third-order polynomial. 

Equation (3) cannot be integrated in closed form; therefore, the Gaussian technique 
of numerical integration was used to evaluate the integral c w dx for the exponential 
curve fit, and the heating rate was computed from equation (1). The heating-rate history 
computed from the integrated exponential expression is presented in figure 6 along with 
the simulated heating rate and the heating rate computed by the method of reference 2 
(third -order -polynomial curve fit). Agreement between the heating rate computed from 
the integrated exponential expression and the simulated heating rate is found to be within 
2 percent at the highest heating rate. 

Although the exponential expression produces more accurate representation of 
heating data, it does not correctly represent the variation of c w through the wall. 
Figure 5 shows an inflection in the c w ~ curve when, in fact, the curve should have 
no points of inflection, (very pronounced inflection points can result when experimental 
scatter is assumed in the four simulated values of c w Also, since heat losses from 

the back wall are assumed to be negligible, the gradient of c w with x should be 
zero at the inside surface. To satisfy these requirements, a procedure was developed 
for determining the restricted third-order exponential function with a least-squares fit to 
the four values of c w The restrictions are that its x derivative be zero at the 
inside surface and that there be no change in the sign of its second x derivative through 
the wall (no inflection point). The development of the method (hereinafter called the 
restricted-exponential method) is presented in the appendix. Results obtained with this 
restricted-exponential method were investigated for the more realistic conditions con- 
sidered in the following section. 


Data Scatter and Sensor Spacing 

Data from free-flight experiments inherently possess some degree of scatter as a 
result of random inaccuracies in the telemetry and data readout systems. Also, the 
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location of the thermal sensors in the wall may differ from flight test to flight test. The 
sensitivity of the restricted-exponential method to these factors was investigated for a 
1.524 -cm -thick beryllium wall. 

Effect of data scatter. - Scattered data (simulated flight temperature data) were gen- 
erated by applying random scatter (produced by a random number generator) to the smooth 
temperature histories of figure 1. A standard deviation of 15° R (8.3° K) was used in the 
investigation. The scattered temperature data were fitted with sixth-order polynomials 
as outlined in the previous section and the derivative histories are shown in figure 7. 
Values of c w at each measurement location were determined and are shown by the 
symbols in figure 8 for t = 2, 6, and 7 seconds. Also shown is the actual distribution of 
c w through the wall at these three times, as in figures 4 and 5. The simulated data 
do not lie on the actual distribution curve, because they were derived from a fairing of 
scattered temperature data. The restricted -exponential least-squares fit of the points 
is a better representation of the actual distribution than is the third-order -polynomial 
curve fit. 

Integration of the c w distributions derived from the scattered temperature 
data produced the heating-rate data of figure 9 for a 1.524 -cm -thick beryllium wall. The 
restricted -exponential method improved the accuracy at the higher heating rates; however, 
for the lower heating rates (q ~ 200 to 400 Btu/ft 2 -sec or 226.98 to 453.96 watts/cm 2 ), 
the restricted-exponential method was less accurate than the third-order -polynomial 
method of reference 2. 

The heating rates computed by both methods oscillate about the simulated heating- 
rate history. These oscillations result from oscillations in the history obtained 

from the sixth -order -polynomial fit of the scattered temperature data. Figure 7 shows 
these oscillations in the history obtained from the scattered data. Since the four 

AT 

values of c w ^ 4 - define the curves which are integrated to obtain heating rates at spec- 
ified times, the oscillations in -2A- directly produce the oscillations in the computed 
heating-rate history. Improved techniques to curve -fit the temperature histories may 
eliminate some of the oscillations and the resulting error in the heating data. 

Effect of sensor spacing. - The number and location of thermal sensors in a thick- 
wall calorimeter vary from one flight test to the next, depending on the experimental 
requirements and the number of telemetry channels available. Project Fire (refs. 2 
and 3) utilized four thermocouples equally spaced through the wall. To determine the 
effects of unequal sensor spacing on the accuracy of the restricted-exponential method, 
the 0.600-inch -thick (1.524-cm) beryllium model was analyzed with the two interior mea- 
surements closer to the heated surface. For the analysis, the two interior temperature 
measurements were located 0.100 and 0.300 inch (0.254 and 0.762 cm) from the heated 
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surface. Such an arrangement would more accurately detect the shape of the temperature 
profile through the wall. 

Temperature data at each thermocouple location were obtained for the simulated 
heating history of figure 1. The scatter applied to the temperature at each thermocouple 
was identical with that applied to the temperatures at the equally spaced thermocouples 
in the previous section. The resulting data were fitted with the restricted-exponential 
expression and compared with the actual distribution and the third-order -polynomial 
method of reference 2. Figure 10 presents the c w |p variation through the wall for 
the unequally spaced thermocouple arrangement. For both thermocouple arrangements, 
the restricted-exponential curves define the gradient through the wall more accurately 
than the third-order polynomial (figs. 10 and 8). 

Heating rates are shown in figure 11 for the unequal thermocouple spacing. The 
restricted-exponential method produces heating rates in close agreement with the exact 
heat input at the highest heating rates. At the lower heating rates, the restricted- 
exponential method also produced more accurate results than the third-order polynomial, 
unlike the results of figure 9 for equal sensor spacing. 

CONCLUSIONS 

The integral method of determining heating rates from temperatures measured at 
four locations within a thermally thick wall has been analyzed for rapidly rising heating 
rates and large thermal as well as physical thicknesses. An exponential function with 
two restrictions has been developed to curve -fit the measured values of c w (T)^E (where 
c w denotes specific heat of wall material; T, temperature; and t, time) and describe 
its variation through the wall. The restricted -exponential method has been analyzed 
for a 1.524-cm-thick beryllium wall with smooth and scattered data and for two sensor 
spacing arrangements. The following conclusions are noted: 

1. The restricted-exponential method provides a better representation of the c w ^ 
distribution through the wall than the previously used polynomial curve fit. 

2. With data scatter (standard deviation of 15° R, or 8.3° K) and equal sensor 
spacing, the restricted exponential method is more accurate than the polynomial method 
at higher heating rates but less accurate at lower heating rates. 
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3. The restricted-exponential method is more accurate than the polynomial method 
for both high and low heating rates for scattered data and unequal sensor spacing (interior 
sensors installed close to the heated surface to detect the slope of the temperature pro- 
file through the wall more accurately). 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., April 24, 1969, 

711-02-09-01-23. 
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APPENDIX 


OPTIMAL FIT OF THIRD-ORDER EXPONENTIAL FUNCTION TO DATA 


A method of obtaining an optimal fit (in the sense of least squares) of the third- 

order exponential c w = e a 0 +a l x+a 2 x +a 3 x ? subject to constraints, to a set of data 
is presented. In order to satisfy the conditions that the c w variation with x should 
have no inflection points and have a zero slope at the inside wail, the following constraints 
have been imposed on the exponential curve fit: 


= 0 at z = 0 


(Al) 


2 

= 0 for all z on [o,£J 
3z 2 


(A2) 


where 


z = l - x 


(A3) 


and 


y = log e c w = a 0 + ajz + a 2 z 2 + a 3 z 3 (A4) 

Expression ( A2) is used in place of -2— (c w ) £ 0, which is the direct requirement that 

dz 2 \ at / 

there be no inflection points in the c w ^ curve, in order to make the procedure more 
tractable. Neglecting gradients of c w yields 




2 


The last term on the right-hand side is positive because it is squared, and ^ is posi- 
tive so that the constraint £ 0 requires that -^r(c w £ 0. It is noted that the 

3z 2 az 2 \ at / 

solution developed herein for an optimal third -order -exponential curve fit will not repre- 
sent the best fit for a quartic exponential. 
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APPENDIX 

The condition of equation (Al) is easily satisfied by setting 

a-i = 0 

which is the reason for introducing z by equation (A3). Then from equation (A4), 

3^v 

— x = 2a 2 + 6a 3 z 
9z z 


(A5) 


Thus, the inequality (A2) may be rewritten as 


a 2 = 0 

2a 2 + 6a 3 £ = 0 


(A6a) 

(A6b) 


In the usual manner of the method of least squares, the sum of the squares of the resid- 
uals S is formed (ref. 8): 


n 2 

S = ^ (a 0 + a 2 z t 2 + a 3Zi 3 - yj (A7) 

i=l 

where z^ is the location of the temperature measurements and y^ is the measured 
value at the corresponding ith locations. 

The problem is reduced to finding a set of values for ag, a 2 , and a 3 that satisfy 
the inequalities (A6) and minimizing S. The sum of the squares of the residuals S may 
easily be written as a quadratic in ag, a 2 , and a 3 , so that the problem is one of qua- 
dratic programing (ref. 9). However, the present problem is sufficiently simple to obviate 
the elegant methods of quadratic programing. 

The first step of the procedure is to minimize S with respect to ag, &2’ anc * a 3* 
Consequently, 

-£S_ _ _9S_ _ _8S_ _ g ( A8 ) 

aaQ 8a2 9^3 

Generating the partial derivatives of equation (A8) from equation (A7) and setting them 
equal to 0 produces 
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na ° + (Ji 


Zl2 ) a2+ (1 



(A9a) 



(A9b) 



H n 

Y 

L=1 



a 2 + 



z i®J a 3 = y i z i 3 


(A9c) 


Equations (A9a) to (A9c) can be easily solved for ag, ag, and ag by the use of 
Cramer's rule. If the resulting values of ag, ag, anc * a 3 satisfy the inequalities (A6), 
the problem is solved. 

If both inequalities (A6a) and (A6b) are violated, the data indicate a cubic for which 

— ^ < 0 over the entire interval. The best fit to these data which satisfies the inequal - 

dz 2 

ities (A6a) and (A6b) is simply 



a 2 = a 3 = 0 


However, it is unrealistic to attempt to fit a curve with positive curvature to data which 
are best fit by a curve with negative curvature throughout the region of interest, and this 
case is rejected. 

In the case of ag < 0 but 2a2 + 6ag£ = 0, it is necessary to set ag = 0 and min- 
imize S in equation (A7) with respect to ag and ag. This leads to equations (A9a) 
and (A9c) with ag = 0. These two linear equations in two unknowns are then solved to 
give ag and ag. If now ag > 0, these values of ag, a -2 (=0), and ag constitute the 
solution. If ag ^ 0, the case is rejected, since the data are again best fit by a curve with 
negative curvature throughout the region of interest. 

Likewise, in the case of a .2 = 0 but 2a2 + 6agZ < 0, it is necessary to set 
2ag + 6ag l = 0, which reduces to 
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APPENDIX 

a 3 = -ff (A10) 

Then, S is minimized with respect to ag and a2, whereby the equation 


3 S 3 S q 

9a 0 9a 3 ~ 


gives 


na o + a 2 


i=i 



yi 


(Alla) 


a 0 


n 


1 

i=l 





(A lib) 


Equations (All) may be solved for ag and a2> and ag is then computed from equa- 
tion (A10). If a2 > 0, the values of a 0 , a. 2 , and a .3 are the acceptable solutions. 

However, if a .2 = 0, the case is rejected as before. 


A schematic diagram of the procedure of fitting the third-order exponential to the 
data is presented in figure 12. 


Proof that the foregoing procedure gives the best fit of the third-order exponential 
to the data under the given restraints can be given geometrically. As was previously 
noted, S is a quadratic in ag, a 2 , and a3. In the ag, a2, a3 space, S = Constant 
describes an ellipsoid which degenerates to a point for ag, a. 2 , and a3 equal to the 
best -fit values computed from equation (A8). At this point, S = S m i n . For any Si and 
S2 with Si > S2 > S m i n , the ellipsoid S = S2 is entirely contained within the ellipsoid 
S = Si, and any point on S2 corresponds to a better fit to the data than any point on Si. 

Figure 13 shows the geometry in the a2~ag plane. For clarity the ag-dimension 
has been suppressed and the ellipsoids are shown as ellipses, but the discussion applies 
to the three-dimensional case. The planes a2 = 0 and 2a2 + 6a 3^ = 0 divide the space 
into four regions. 

Region I is denoted as the set of points ^ag, a2, a 3) for which the inequalities (A6) 
are satisfied, as indicated by point A. Region n is the set of points for which neither of 
the inequalities (A6) are satisfied, as shown by point B. Region m is the set of points 
for which 
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a 2 < ° 


2a£ + 6agZ 1= 0 


as given by example points C, D, and E. Region IV is the set of points for which 

a 2 =0 

2a2 + 6agZ < 0 


as given by example points F and G. 

If the solution of equations (A9) is in region I, it satisfies all requirements and is 
accepted. If the solution is in region n, the data are best fit by a cubic with completely 
negative curvature and the case is rejected because a least-squares fit with the con- 
straints of positive curvature would be meaningless in application. If for the given data 
the solution of equations (A9) falls in region HI, then equations (A9) do not apply. In this 
case, the (ag, a 2 , a 3 ^ point for which S is a minimum and which satisfies the inequal- 
ity (A6a) is the point of tangency of an ellipsoid which is tangent to the plane a 2 = 0, such 
as point in figure 13. This point is given by the solution to equations (A9a) and (A9c) 
with a 2 = 0. However, this a 2 = 0 point may fall on the boundary of region II, such as 
point Dj, where ag ^ 0. In this event, the case is rejected. 

Finally, as with point E, the smallest ellipsoid intersecting region I (touching 
region I at point Ej) may also intersect region n, so that although the given procedure 
gives a best fit satisfying all conditions (positive curvature throughout the interval), a 
better fit may be found with negative curvature throughout the interval. If the solution 
of equations (A9) lies in region IV, the same type of considerations apply as when the 
solution lies in region III. 
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Figure 1.- Heating-rate history and associated temperature history at different locations in 1.524-cm-thick beryllium wall. 


Temperature 




q, Btu/ft -sec 



Time, sec 

Figure 3.- Simulated heating rate and heating rate computed with third-order polynomial curve fit. 
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Figure 6.- Simulated heating rate and heating rates computed by integral method with third-order polynomial and 

third-order-exponential curve fits. 
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Figure 7.- Temperature-derivative histories (from sixth-order polynomial) for smooth and scattered data. 
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Figure 11.- Heating-rate histories for unequal sensor spacing. 
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